Triangulating stars on the night sky

In a previous notebook, I have shown properties of the distribution of stars in the sky. Here, I would like to use the existing database of stars' positions and display them as a triangulation

Let's first initialize the notebook:

In [1]:
import numpy as np
np.set_printoptions(precision=6, suppress=True)
%matplotlib inline
%config InlineBackend.figure_format = 'svg'
import matplotlib.pyplot as plt
phi = (np.sqrt(5)+1)/2
fig_width = 10
figsize = (fig_width, fig_width/phi)

importing data

Importing all stars' position is as simple as invoking the HYG database:

In [2]:
import pandas as pd
url = "https://github.com/astronexus/HYG-Database/raw/master/hygdata_v3.csv"
space = pd.read_csv(url, index_col=0)
In [3]:
print(f'Columns in the database: {space.columns=}')
Columns in the database: space.columns=Index(['hip', 'hd', 'hr', 'gl', 'bf', 'proper', 'ra', 'dec', 'dist', 'pmra',
       'pmdec', 'rv', 'mag', 'absmag', 'spect', 'ci', 'x', 'y', 'z', 'vx',
       'vy', 'vz', 'rarad', 'decrad', 'pmrarad', 'pmdecrad', 'bayer', 'flam',
       'con', 'comp', 'comp_primary', 'base', 'lum', 'var', 'var_min',
       'var_max'],
      dtype='object')
In [4]:
print(f'Number of stars in the catalog = {len(space)=}')
Number of stars in the catalog = len(space)=119614

extracting ra, dec and mag

For which we may extract what interests us: position (right ascension and declination) and visual magnitude:

In [5]:
space_pos = space[['ra', 'dec', 'mag']]
space_pos
Out[5]:
ra dec mag
id
0 0.000000 0.000000 -26.70
1 0.000060 1.089009 9.10
2 0.000283 -19.498840 9.27
3 0.000335 38.859279 6.61
4 0.000569 -51.893546 8.06
... ... ... ...
119611 23.963895 38.629391 12.64
119612 23.996567 47.762093 16.10
119613 23.996218 -44.067905 12.82
119614 23.997386 -34.111986 12.80
119615 0.036059 -43.165974 13.05

119614 rows × 3 columns

First, right ascension is "the angular distance of a particular point measured eastward along the celestial equator from the Sun at the March equinox to the (hour circle of the) point in question above the earth." It is given in hours:

In [6]:
ra_min, ra_max = 0, 24
print(f"RA: {space_pos['ra'].min()=}, {space_pos['ra'].max()=}")
RA: space_pos['ra'].min()=0.0, space_pos['ra'].max()=23.998594

normalizing ra into az

Let's convert this in visual angle, that is in an azimuth:

In [7]:
space_norm = space_pos[['dec', 'mag']].copy()
space_norm
Out[7]:
dec mag
id
0 0.000000 -26.70
1 1.089009 9.10
2 -19.498840 9.27
3 38.859279 6.61
4 -51.893546 8.06
... ... ...
119611 38.629391 12.64
119612 47.762093 16.10
119613 -44.067905 12.82
119614 -34.111986 12.80
119615 -43.165974 13.05

119614 rows × 2 columns

In [8]:
az_min, az_max = 0, 360
def ra2az(ra):
    return az_max - ra / ra_max * az_max
In [9]:
space_norm["az"] =  ra2az(space_pos['ra'])
In [10]:
print(f"AZ: {space_norm['az'].min()=}, {space_norm['az'].max()=}")
AZ: space_norm['az'].min()=0.021090000000015152, space_norm['az'].max()=360.0

Then, declination is "comparable to geographic latitude, projected onto the celestial sphere" and is given here in degrees:

In [11]:
dec_min, dec_max = -90, 90
print(f"DEC: {space_norm['dec'].min()=}, {space_norm['dec'].max()=}")
DEC: space_norm['dec'].min()=-89.782428, space_norm['dec'].max()=89.569427

The magnitude varies a lot (the less the value, the more it is visible - "The apparent magnitudes of known objects range from the Sun at −26.7 to objects in deep Hubble Space Telescope images of magnitude +31.5"):

In [12]:
print(f"mag: {space_norm['mag'].min()=}, {space_norm['mag'].max()=}")
mag: space_norm['mag'].min()=-26.7, space_norm['mag'].max()=21.0

Let's normalize the lower bound:

In [13]:
space_norm['mag'] = space_pos['mag'] - space_pos['mag'].min()
print(f"mag: {space_norm['mag'].min()=}, {space_norm['mag'].max()=}")
mag: space_norm['mag'].min()=0.0, space_norm['mag'].max()=47.7

Let's define a threshold using the quantile function:

In [14]:
print(f"{space_norm['mag'].quantile(q=0.01)=}")
space_norm['mag'].quantile(q=0.01)=31.45

Which leaves only a limited number of stars

In [15]:
space_bright = space_norm[space_norm['mag']<space_norm['mag'].quantile(q=0.05)]
In [16]:
#space_bright = space_pos[space_pos['mag']<6]
In [17]:
print(f'Number of bright stars = {len(space_bright)=}')
Number of bright stars = len(space_bright)=5955
In [18]:
print(f"MAG: {space_bright['mag'].min()=}, {space_bright['mag'].max()=}")
MAG: space_bright['mag'].min()=0.0, space_bright['mag'].max()=32.85

scatter plots

From these elements, we may plot the stars on these coordinates:

In [19]:
stars_color = np.array([255., 235., 220.])
stars_color /= 255.
In [35]:
fig, ax = plt.subplots(figsize=(fig_width, fig_width/phi))
ax.set_facecolor('black')
ax.scatter(space_bright['az'], space_bright['dec'], c=stars_color[None, :], ec='none', s=.1 * (space_bright['mag'].max()-space_bright['mag']))
ax.set_xlabel('Azimuth')
ax.set_ylabel('Declination')
ax.set_xlim(az_min, az_max)
ax.set_ylim(dec_min, dec_max);
2021-12-04T18:17:27.159042image/svg+xmlMatplotlib v3.4.3, https://matplotlib.org/

do you see the Milky way?

Orion

We may want to focus on the Orion constellation):

In [21]:
orion_az_max, orion_az_min = ra2az(4.8), ra2az(6.25)
orion_dec_min, orion_dec_max = -13.0, 22.
space_orion = space_bright[(orion_az_min < space_bright['az']) * (space_bright['az'] < orion_az_max) *
                                (orion_dec_min < space_bright['dec']) * (space_bright['dec'] < orion_dec_max)]
In [22]:
print(f'Number of bright stars in orion = {len(space_orion)=}')
Number of bright stars in orion = len(space_orion)=192
In [23]:
print(f"MAG: {space_orion['mag'].min()=}, {space_orion['mag'].max()=}")
MAG: space_orion['mag'].min()=26.88, space_orion['mag'].max()=32.85
In [24]:
fig, ax = plt.subplots(figsize=(fig_width, fig_width))
ax.set_facecolor('black')
ax.scatter(space_orion['az'], space_orion['dec'], c=stars_color[None, :], s= 10 * (space_orion['mag'].max()-space_orion['mag']))
ax.set_aspect('equal')
ax.set_xlabel('Azimuth')
ax.set_ylabel('Declination')
ax.set_xlim(orion_az_min, orion_az_max)
ax.set_ylim(orion_dec_min, orion_dec_max);
2021-12-04T17:57:17.162736image/svg+xmlMatplotlib v3.4.3, https://matplotlib.org/

Which is to be compared to this image:

In [25]:
from  IPython.display import display, SVG, Image
display(SVG('https://upload.wikimedia.org/wikipedia/commons/f/ff/Orion_IAU.svg'))

Triangulation

Following work with Etienne Rey, we have developped visualization based on triangulation of an optimal packing of point clouds.

In [26]:
display(Image('https://laurentperrinet.github.io/post/2021-10-04_interstices/featured.jpg'))

The idea here is to apply such a visualization to the stars position. IT is based on

In [27]:
import matplotlib.tri as tri
triang = tri.Triangulation(space_bright['az'], space_bright['dec'])
In [28]:
fig, ax = plt.subplots(figsize=(fig_width, fig_width))
ax.set_facecolor('black')
ax.triplot(triang, lw=0.1, color=stars_color[None, :])
ax.set_aspect('equal')
ax.set_xlabel('Azimuth')
ax.set_ylabel('Declination')
border_ratio = .95
ax.set_xlim(az_min+(az_max-az_min)*(1-border_ratio), az_min+(az_max-az_min)*border_ratio)
ax.set_ylim(dec_min+(dec_max-dec_min)*(1-border_ratio), dec_min+(dec_max-dec_min)*border_ratio);
2021-12-04T17:57:23.664950image/svg+xmlMatplotlib v3.4.3, https://matplotlib.org/
In [29]:
fig, ax = plt.subplots(figsize=(fig_width, fig_width))
ax.set_facecolor('black')
ax.triplot(triang, lw=0.3, c=stars_color)
ax.scatter(space_orion['az'], space_orion['dec'], c=stars_color[None, :], ec='none', alpha=.5, s= 40 * (space_orion['mag'].max()-space_orion['mag']))
ax.set_aspect('equal')
ax.set_xlabel('Azimuth')
ax.set_ylabel('Declination')
ax.set_xlim(orion_az_min+(orion_az_max-orion_az_min)*(1-border_ratio), orion_az_min+(orion_az_max-orion_az_min)*border_ratio)
ax.set_ylim(orion_dec_min+(orion_dec_max-orion_dec_min)*(1-border_ratio), orion_dec_min+(orion_dec_max-orion_dec_min)*border_ratio);
2021-12-04T17:57:24.087940image/svg+xmlMatplotlib v3.4.3, https://matplotlib.org/
In [30]:
q = 0.01
space_norm['mag'].quantile(q=q)
Out[30]:
31.45
In [31]:
space_orion = space_norm[space_norm['mag']<space_norm['mag'].quantile(q=q)]
In [32]:
space_orion = space_orion[(orion_az_min < space_orion['az']) * (space_orion['az'] < orion_az_max) *
                           (orion_dec_min < space_orion['dec']) * (space_orion['dec'] < orion_dec_max)]

Let's now super-impose different triangulations at different magnitudes:

In [33]:
fig, ax = plt.subplots(figsize=(fig_width, fig_width))
ax.set_facecolor('black')

for mag in range(2, 8):
    #q = 1 - 0.8**mag
    #print (q)
    space_orion = space_norm.copy()
    #space_orion = space_norm[space_norm['mag']<mag]
    space_orion = space_orion[(orion_az_min < space_orion['az']) * (space_orion['az'] < orion_az_max) *
                               (orion_dec_min < space_orion['dec']) * (space_orion['dec'] < orion_dec_max)]
    space_orion = space_orion[space_orion['mag'] > space_orion['mag'].max()-mag]
    print(f'{mag=:.3f}, Number of bright stars = {len(space_orion)=}')
    triang = tri.Triangulation(space_orion['az'], space_orion['dec'])
    ax.triplot(triang, lw=9/mag**2, color=stars_color, alpha=.4)
#ax.scatter(space_orion['az'], space_orion['dec'], s= 40 * (space_orion['mag'].max()-space_orion['mag']))
ax.set_aspect('equal')
ax.set_xlabel('Azimuth')
ax.set_ylabel('Declination')
ax.set_xlim(orion_az_min+(orion_az_max-orion_az_min)*(1-border_ratio), orion_az_min+(orion_az_max-orion_az_min)*border_ratio)
ax.set_ylim(orion_dec_min+(orion_dec_max-orion_dec_min)*(1-border_ratio), orion_dec_min+(orion_dec_max-orion_dec_min)*border_ratio);
mag=2.000, Number of bright stars = len(space_orion)=8
mag=3.000, Number of bright stars = len(space_orion)=14
mag=4.000, Number of bright stars = len(space_orion)=25
mag=5.000, Number of bright stars = len(space_orion)=51
mag=6.000, Number of bright stars = len(space_orion)=140
mag=7.000, Number of bright stars = len(space_orion)=452
2021-12-04T17:57:24.657106image/svg+xmlMatplotlib v3.4.3, https://matplotlib.org/

some book keeping for the notebook

In [34]:
%load_ext watermark
%watermark -i -h -m -v -p numpy,matplotlib,scipy,pillow,imageio  -r -g -b
Python implementation: CPython
Python version       : 3.9.9
IPython version      : 7.30.0

numpy     : 1.20.3
matplotlib: 3.4.3
scipy     : 1.6.0
pillow    : not installed
imageio   : 2.13.1

Compiler    : Clang 13.0.0 (clang-1300.0.29.3)
OS          : Darwin
Release     : 21.1.0
Machine     : x86_64
Processor   : i386
CPU cores   : 4
Architecture: 64bit

Hostname: ekla

Git hash: 6b0c86af026462cfae61f02aef25a3acb568faf9

Git repo: https://github.com/laurentperrinet/sciblog.git

Git branch: master

Contents © 2022 Laurent Perrinet - Powered by Nikola.
Creative Commons License BY-NC-SA